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ABSTRACT 

It has been recently found that the local fluctuations of the QSO's Lya absorption 
spectrum transmitted flux show spiky structures. This implies that the mass fields of 
the intergalactic medium (IGM) is intermittent. This feature cannot be explained by 
the clustering evolution of cosmic mass field in the linear regimes and is also difficult 
to incorporate into the hierarchical clustering scenario. We calculate the structure 
functions and intermittent exponent of the IGM and HI for full hydrodynamical 
simulation samples. The result shows the intermittent features of the Lya transmitted 
fiux fluctuations as well as the mass field of the IGM. We find that within the error 
bars of current data, all the intermittent behavior of the simulation samples are 
consistent with the observation. This result is different from our earlier result (Pando 
et al 2002), which shows that the intermittent behavior of samples generated by 
pseudo-hydro simulation cannot be fltted with observed data. One difference between 
the pseudo-hydro and full hydro simulations is in treating the dynamical relation 
between the IGM (or HI) and dark matter fields. The former assumes that the IGM 
density distribution traces the underlying dark matter point-by-point on scales larger 
than the Jeans length in either the linear or nonlinear regimes. However, hydrodynamic 
studies have found that a statistical discrepancy between the IGM field and underlying 
dark matter in nonlinear regime is possible. We find that the point-by-point correlation 
between the IGM density perturbations and dark matter become weaker on comoving 
scales less than 2 h~^ Mpc (in LCDM model), which is larger than the IGM Jeans 
length. 



Subject headings: cosmology: theory - large-scale structure of the universe 

^Center for Astrophysics, University of Science and Technology of China, Hefei, Anhui 230026, P.R.China 

^National Astronomical Observatories, Chinese Academy of Science, Chao-Yang District, Beijing, 100012, P.R. 
China 

^Department of Physics, DePaul University, Chicago II, 60614 
^Department of Physics, University of Arizona, Tucson, AZ 85721 



- 2 - 



1. Introduction 

In the past few years there has been tremendous progress in understanding the evolution 
of cosmic structures on scales from galaxies to about 100 h~^ Mpc. This progress has occurred 
mainly in two distinct areas: 1.) constraining cosmological parameters such as the mass density, 
dark energy, power spectrum of the initial perturbations of the cosmic mass field, etc., 2.) the 
modeling of the origin of light-emitting objects to explain the observed mass density profiles, 
velocity dispersions, mass-to-light ratios, substructures, bias behavior, etc., of quasars, galaxies 
and galaxy clusters. 

These two distinct areas coincide with two different regimes in the context of the dynamics 
of structure formation. The first one is based mainly on the evolution of the cosmic mass field 
in the linear approximation. The basic statistical tool is the power spectrum, either observed or 
recovered by Gaussianization, of the weakly non-linear field. The second regime is characterized 
by the dynamics of collapsed massive halos. The formation and evolution of luminous clumps can 
be modeled by the growing and merging of massive halos from initially Gaussian perturbations of 
the mass and velocity fields. The hierarchical clustering of halos is the most effective scenario for 
modeling galaxy formation. 

However, the clustering behavior of the mass and velocity fields of Inter-Galactic Medium 
(IGM), i.e., the major component of baryonic matter in the universe, is missing from these two 
regimes. It is certain that the IGM field is not linear and Gaussian because the distributions of 
QSO's Lya forests from either observational data or simulation samples (Bi 1993, Bi, Ge & Fang 

1995) , have been found to be non-Gaussian even on scales as large as a few h~^ Mpc (Pando &; 
Fang 1996, 1998a). Neither is the evolution of the IGM field well modeled by the dynamics of 
collapsing halos, as the QSO's Lya forests are generally not gravitationally confined (Fang et al 

1996) . The success of modeling the QSO's Lya forests via a lognormal model of the IGM (Bi & 
Davidsen 1997) directly indicates that the IGM field probably does not (or has not yet) undergone 
a correlation hierarchy evolution. 

The clustering behavior of the IGM mass field can play an important role in the dynamical 
gap between linear power spectrum and collapsed massive halos. For instance, it has recently 
been found that the QSO Lya transmitted flux fluctuations are intermittent (Jamkhedkar, Zhan 
&; Fang 2000, Feng, Pando &; Fang 2001, Zhan, Jamkhedkar & Fang 2001). A field is intermittent 
if the 2n*'* power of the local differences grows faster than the n*'' power of the variance of these 
differences (see Section 2.) This result cannot be found using the statistical tools of the linear 
regime such as power spectrum. At the same time, it is unclear how to incorporate intermittent 
behavior with the hierarchical clustering scenario. The intermittent features have also been studied 
with simulated Lya forests produced by a pseudo hydro technique. All the pseudo hydro samples 
for low density cold dark matter (LCDM) and warm dark matter (WDM) models are found to 
be substantially intermittent. Although the intermittent features of the simulation samples were 
different from those of the observations (Pando et al 2002), the results implies that intermittency 
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is a common feature of non-linear evolution and is potentially useful for model discrimination. 

This motivated us to re-do the intermittent analysis with simulated Lja transmitted fluxes 
produced by a full hydrodynamical treatment. This motivation was strengthened when we realized 
that the relation between the IGM and dark matter fields in the non-linear regime is treated 
differently in the pseudo-hydro and full hydro simulations. In the pseudo-hydro simulation, the 
density of the IGM in each pixel is assigned according to the dark matter density at that pixel. 
That is, the pseudo-hydro simulation assumes that the density distribution of IGM traces the 
underlying dark matter point-by-point on scales larger than the Jeans length in either the linear 
or nonlinear regimes. This assumption looks reasonable when one considers that the IGM is only 
a passive substance in the cosmic mass field since its evolution is dominated by the gravity of dark 
matter. 

However, hydrodynamic studies have shown that a passive substance in the system docs not 
always behave like the underlying dominant mass field when non-linear evolution has taken place 
(for a review, Shraiman &; Siggia 2001). The statistical properties of a passive substance can 
decouple from those of the underlying field. As we show here, the dynamical equations of the IGM 
within the gravitational field of dark matter are the same as that used for describing the statistical 
decoupling between the passive substance and underlying field. Therefore, it is possible there is a 
statistical discrepancy between the IGM field and underlying dark matter in nonlinear regime. 

The purpose of this paper is then two-fold. First we study whether the intermittency of the 
QSO Lya transmitted flux is matched by the full hydro simulation samples. Secondly we search 
for possible signs of the statistical discrepancy between the IGM field and the underlying dark 
matter with these samples. The paper is organized as follows. §2 presents a brief introduction 
to intermittent fields and the method to quantify intermittent random fields. §3 describes the 
observed samples that we use. §4 describes full hydrodynamic simulations. The statistical 
deviation of the IGM distribution from the dark matter field is shown in §5. A comparison between 
the model predicted and observed intermittency is performed in §6. Finally, the conclusions and 
discussions are presented in §7. The dynamical equations for the IGM evolution in non-linear 
regime is given in Appendix. 

2. Intermittency 

Much of this section has been presented in our previous work (Pando et al 2002). For the 
paper to be self-contained, we repeat these results briefly. 
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2.1. The intermittent exponent 

Intermittency is used to characterize a random field with 'spike-gap-spike' features, i.e., 
structures that are essentially strong enhancements randomly and widely scattered in space and/or 
time, with a low field value between the spikes. The spike-gap-spike feature is more pronounced 
on smaller scales. This feature was originally foTind in the temperature and velocity distributions 
in turbulence (Batchelor & Townsend 1949; Prisch 1995). 

An intermittent random density field /9(x) is defined by the asymptotic behavior (indicated 
below by x) of the ratio between the high- and low-order moments of the field as (Gartner &; 
Molchanov, 1990; Zel'dovich, Ruzmaikin, & Sokoloff, 1990) 



([/9(x + r)-Kx)] 



2n\ / ^ \ C 



(1) 



[([Kx + r)-p(x)]2)]n \L 

where L is the size of the sample. Intermittency is measured by the n- and r-dependcncics of (" - 
the intermittent exponent. If the exponent C is negative on small scales r, the field is intermittent. 
That is, the ratio in eq. (1) for an intermittent field diverges as r ^ 0. 

This divergence cannot be measured by individual higher order correlation functions. 
Although three- and four-point correlation functions are useful for distinguishing a Gaussian from 
a non-Gaussian field, they are insensitive to the difference between an intermittent field and 
a non-Gaussian, but non-intermittent field. To measure intermittency, we define the structure 
function as 

S^- = {\Ar{x)n (2) 

where Ar(x) = p{x + r) — p{x), and (...) is the average over the ensemble of fields. If the field 
is statistically homogeneous, S^.^ is independent of x and depends only on r. 5^" is called the 
structure function. When n = 1, we have 

S', = {\Ar{x)f). (3) 

is the mean of the square of the density fluctuations at wavenumber k ~ 27r/r, and therefore, 
is the power spectrum of the field. With the structure function, eq.(l) can be rewritten as 



Q2n 



r 



r 



The ratio S^."' / {S^.Y'' is the n*^ moment S*^" normalized by the power S'^. Since measures the 
width of the probability distribution function (PDF) of |Ar(x)|, while S"^" is sensitive to the tail 
of the PDF, the ratio eq.(4) measures the fraction of density fluctuations in the tail on different 
scales. 

For an intermittent field the ratio of S^"' to [S^]"' is larger for smaller r and indicates that the 
field contains "abnormal" events of large density fluctuations |Ar(a;)|. That is, the PDF of |Ar(x)| 
is long-tailed. The long-tail events correspond to a sharp increase or decrease of the density /9(x). 
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They are rare events and yield high peaks or spikes in the distribution of density difference. 

Generally speaking, the intermittent exponent ^ measures the smoothness of the field: for positive 
the field is smoother on smaller scales. If is negative, the field is rough on small scales. 

The basic quantity used for calculating the intermittent exponent ( is the density difference 
|Af.(x)|. The exponent ( docs not rely on the mean density p and can be measured even if the 
mass field is fractal, in which case the mean density is no longer available. 



2.2. Examples of the intermittent exponent 

In our previous work we presented several examples of intermittency in various kinds of fields. 
Here we summarize those results for the kinds of fields that pertain to this work. 

For a Gaussian field, 

q2n 

This ratio is independent of scale r, and therefore, the intermittent exponent ^ = 0. A Gaussian 
field is not intermittent. A self-similar field is also not intermittent. 

For a hierarchically clustered field, it has been shown (Feng, Pando, &; Fang 2001) that 

g2n . . -(d-«)(n-l) 

oc(-l , (6) 



where d is the spatial dimension, and the coefficient At is a constant depending on the power law 
index of the power spectrum. In the case where k < d, the field is intermittent with exponent 

C~-(d-«;)(n-l). (7) 
Finally, for a lognormal field the PDF of Aj.(x) is given by 



23/27rV2|A^(x)|a(r) { 2 V (^(r) 

where |Ar(x)|Tn, is the median of |Ar(x)| (Vanmarcke 1983). The variance a{r) of ln|Ar(x)| 
generally can be a function of the scale r. Using eq.(8), we have 

[52]- - e • lyj 

The intermittent exponent of a lognormal field is then 

C ~ 2(n2 -n)c72(r)/ln(r/L). (10) 
Because r < L, ^ is negative. Therefore, a lognormal field is intermittent. 
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To summarize, the structure function and intermittent exponent provide a complete and 

unified description of intermittent fields. The n- and r-dependencies of the structure functions 
and intermittent exponent (" are sensitive to the details of the intermittency of the field. These 
measures are very powerful for distinguishing among fields that are Gaussian, self-similar, 
mono-fractal, multi-fractal, and singular. 



2.3. The intermittent exponent in the wavelet basis 

The quantity Ar{x) or [p{x + r) — p{x)] contains two variables: the position x and the scale 
r, and therefore, p{x + r) — p{x), can best be calculated with a space-scale decomposition. Among 
the mathematical tools used to decompose functions in scale space as well as physical space, the 
discrete wavelet transform (DWT) is an excellent choice. The DWT basis are orthogonal and 
complete, and therefore, the DWT decomposition does not lose information (completeness), and 
does not cause false correlations due to redundancy (non-orthogonality). These two points are 
essential in studying the higher order statistical features of random fields. 

We restrict our discussion to a 1-D random field p{x) extending in a physical or redshift 
spaces L = X2 — xi. To apply the DWT, we first chop the spatial range L = X2 — xi into 2^ 
subintervals labeled I = 0, ...2^ — 1 where j is a positive integer. Each subinterval spans a spatial 
range with the Z*'* subinterval spanning from xi + to xi + L{1 + l)/2-'. This operation 

decomposes the space L into cells (j, /), where j denotes the scale L/2-', and I the spatial range. 

Corresponding to each cell, there is a scaling function (t)ji{x), and a wavelet function ipji{x). 
The scaling function (j)j^i{x) is a window function on scale j and at position I. The wavelet 
functions ^ji{x) form the basis for the scale-space decomposition. The most important property 
of the DWT basis is its locality in both scale and physical spaces. 

With the DWT, a density field p{x) can be decomposed as (Fang &; Thews 1998). 

1=0 j'=j 1=0 

where J is given by the finest scale (resolution) of the sample, i.e., Ax = , and j is the scale 
of interest. The scaling function coefficient (SFC) eji in eq.(ll) is given by projecting p{x) onto 

= J pix)(pj^i{x)dx. (12) 
The SFC Eji describes the mean field of the mode {j, I). 

The wavelet function coefficient (WFC), e^y, in eq.(ll) is obtained by projecting p{x) onto 

= / p{x)tpj,i{^)dx. (13) 
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The wavelet function tpji{x) is admissible (Daubechies 1992), i.e., / il)ji{x)dx = 0, so the WFC is 
basically the difference between p{x) and p{x + r), where x ~ xi + and r L/2K Therefore, 

the WFC eq.(13) measures the fluctuation on scale j and at position I. Because the admissibility 
of wavelets, eq.(13) can also be written as 



j 5{x)il)j^i{x)dx, (14) 



if the normalization is p = 1. Thus, the WFC can be used as the variable p{x + r) — p{x) or 
8{x + r) — 5{x) in eq.(2). All compactly supported wavelet bases produce similar results. We will 
use Daubechies 4 (Daubechies, 1992) in the study below. 

When the "fair sample hypothesis" is applicable (Peebles 1980), the structure function eq.(2) 
can be calculated as the spatial average 

where L is the spatial range of the sample. In the DWT representation, eq.(15) yields 

where j plays the same role as r in eq.(2). S*?"" is the mean of moment le^yP" over the position 
index I. The intermittent exponent C, in eq.(4) can be calculated from 

Generally, C, depends on n and j. 



3. The Observed Lya transmitted flux 

3.1. QSO's Lya absorption spectra 

The observed transmitted flux at wavelength A is F{X) = F^e^'^^-^^ where t(A) is optical 
depth, and Fc the continuum which is generally A-dependent. F^, is given by the emission spectrum 
of the QSO considered. To study the intermittent behavior of the IGM we focus on the fluctuation 
field e""^*^^). The intermittent features of transmitted flux of QSO's Lya absorption spectrum were 
previously studied with a set of about 30 Keck HIRES spectra (Jamkhedkar, Zhan, &; Fang, 2000, 
Pando et al. 2002). This work will consider the two best absorption spectra in that the sample set. 

1. HS1700+6416 

The absorption spectrum of QSO HS1700+6416 (redshift z = 2.72) covers a wavelength range 
from 3727.012 A to 5523.554 A, for a total of 55882 pixels. On average, a pixel is about 0.028 A. For 
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all pixels in this data set, the ratio A A/ A ~ 7 x 10~^ and constant, or in terms of the local velocity, 
5v = cAA/A ~ 2.1 km/s, and therefore, the resolution is about 4 km/s. The distance between 
N pixels in units of the local velocity scale is given by At; = 2c(l — exp [— (1/2)A'^ 6v/c]) km/s, 
or wavenumber k = 27r/Au s/km. As 6v/c <C 1, we have Av ~ N Sv. In this work, we use the 
data from A =3815.6 A to 4434.3 A, which corresponds to z = 2.14 ~ 2.65. The lower limit of 
the wavelength is set to exclude Ly/3 absorption, the upper to eliminate proximity effects. The 
continuum Fc for this spectrum was done by IRAF CONTINUUM fitting by Kirkman and Tytler 
(1997). 

In each pixel the data contain wavelength Aj, flux -F(Aj) and noise a{Xi), which accounts 
for the Poisson fluctuations in the photon count, the noise due to the background and the 
instrumentation. The S/N ratio is about 8. 

This data set contains a fitted flux which smooths out all fluctuations on smallest scales, i.e. 
one or two pixels. The fluctuations on these scales are substantially suppressed by the fitting 
procedure. However, we analyze only scales greater than or equal to 16 pixels. Therefore, the 
smoothing has little effect on our results. Metal lines for this spectrum have been identified down 
to the Doppler parameters 6 ~ 10 km/s by Wei Zheng (unpublished). 

2. qOOU+8118 

This spectra has the highest S/N--^ 30 among the 30 Keck QSO spectra collected by 
Kirkman and Tytler (1997). The redshift of q0014+8118 is 3.387. For all pixels in this data set 
AA/A ~ 13.8 X 10~^, or 5v = cAA/A ~ 4.01 km/s, and, therefore, the resolution is about 8 km/s. 
In each pixel, the data also contains wavelength Aj, flux -F(Aj) and noise a{Xi). The continuum of 
the spectrum is given by IRAF CONTINUUM fitting as well. 

We use the wavelength region from the Ly/? emission to the Lya emission, excluding a region 
of about 0.06 in redshift close to the quasar to avoid any proximity effects. This corresponds to 
z = 2.702 to 3.314, which contains 12703 pixels. The distance between N pixels in the units of the 
local velocity scale is given by Av ~ x 4.01 km/s, or wavenumber k = 27r/Au s/km. Absorption 
features with 6 < 20 km s~^ of this spectrum have been identified. 

3.2. Treatment of unwanted data 

To study the intermittency of the transmission flux fluctuations, we should remove the effects 
of unwanted data, such as 1.) bad pixels (gaps without data), 2.) metal lines, and 3.) negative 
flux pixels. The last category generally consists of saturated absorption regions having lower S/N. 
Although the percentage of low S/N data is not large, it will introduce large uncertainties in the 
analysis. 

The conventional technique for removing unwanted data is to eliminate the unwanted bins and 
to smoothly rejoin the rest of the forest spectra. However, by taking advantage of the localization 
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of wavelets we can use the algorithm of DWT de-noising by thresholding (Donoho 1995) as follows 



1. Calculate the SFCs for both the transmission F{X) and noise (t(A), i.e. 



jl = I F{x)(t>ji{x)dx, efi=l a{x)(t)ji{x)dx. 



(18) 



2. Identify unwanted mode (j, I) using the condition 



<f 



(19) 



where / is a constant. This condition flags all modes with S/N less than /. We can also flag 
modes dominated by metal lines. 

3. Since all the statistical quantities in the DWT representation arc based on an average over 
the modes we do not count all the flagged modes when computing these averages. 



This is the conditional-counting method to treat the unwanted data (Pando et al. 2002). 

It should be emphasized that the condition (19) is applied at each scale j. For instance, if the 
size of a saturated absorption region is d, eq.(19) requires the removal of all modes on scales 
less than d and located within this region, as these modes have very low S/N (ej^/ej^ < /) and the 
local fluctuations in these regions are not reliable. However, condition (19) does not require the 
removal of the saturated absorption region when d is less than the scale of modes (j, I) considered. 
If c-fi/^ji > /) the mode is counted in the statistics, regardless of whether the mode 
contains the saturated absorption regions. The DWT analysis allows an effective application of 
data containing low S/N regions. That is, with the scale-by-scale denoising, it eliminates low S/N 
modes on small scales and in these regions, but keeps their contributions to large scale modes. 
With this procedure, no eliminating and smoothly rejoining of the data is needed, and therefore, 
it avoids corrupting long wavelength modes. 

With this method, we can still calculate the structure functions by eq.(17), but the average 
is not over all modes I, but over the un- flagged modes only. We can also flag modes {j, I) one 
by one which are dominated by metal lines. Since the DWT calculation assumes the sample is 
periodized, this may cause uncertainty at the boundary. To reduce this effect, we drop 3 modes 
near the boundary, and also flag two modes around a unwanted mode to reduce any boundary 
effects of the chunks. 



3.3. Power spectrum of HS1700 



To check the quality of the data, we calculate the power spectrum of the transmitted flux 
fluctuations, AF = F{X) - {F{X)), of HS1700-I-6416. In the DWT representation, this power 
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spectrum is given by (Pando & Fang 1998b; Fang & Feng 2000, Jamkhedkar, Bi, &; Fang, 2001) 

^. = ^EW-^EW.)'- (20) 

The first term on the r.h.s. of eq.(20) is the power of AF on scale j, in which the wavelet 
coefficients (WFC) are given by 

efi = J Fix)^ji{x)dx. (21) 
The second term on the r.h.s. of eq.(20) is due to the noise field a{X) and calculated by 

(e7,)2 = J a\x)i^l{x)dx. (22) 

The DWT power spectrum Pj of a random field F{x) is related to its Fourier power spectrum 
P(n) by 

1 oo 

P3 = ^ E \i^{n/2^)?P{n). (23) 

n=— oo 

where ■0 is the Fourier transform of the wavelet. This implies that the DWT power spectrum Pj 
is the banded Fourier power of the flux fluctuations F{\) — {F{\)). Since for D4 wavelet |V'(r/)p is 
peaked at 77 2± 1, the band j corresponds to the wavenumber k = 27rn/L ~ 27r2^ /L. 

Figure 1 plots the results of the power spectra for HS1700+6416, in which the data is treated 
three ways, 1.) using the original fitted flux, 2.) using the conditional counting method eq.(19) 
with / = 1, and 3.) with / = 3. 

At the first glance, the conditional-counting condition eq.(19) would seem to preferentially 
miss-count modes in the low transmission regions and lead to an /-dependence in the power 
spectrum. However Fig. 1 shows that the power spectrum Pj is independent of / on entire the 
scale range considered from / = 1 to 3. The reason can be seen from eq.(20), which shows that 
the contribution to the power Pj given by mode {j, I) is (ej;)^ — (Cj/)^- The noise subtraction term 
{eji)^ guarantees that the contribution of modes with small ratio S/N to Pj is always small or 
negligible. In other words, all the miss— counted modes make very small or negligible contributions 
to Pj regardless the parameter /. 

Also in Fig. 1, on scales j > 10, or Av < 32 km/s, the power of the fitted flux is lower than 

the power spectrum. This discrepancy is probably because the fitted flux is over-suppressed, i.e., 
not only noise, but also real fluctuations are suppressed by the fitting procedure. This discrepancy 
indicates that the effects of noise become significant on scales Av < 32 km/s. For this reason, we 
use only data on scales Av > 32 km/s. Sample q0014-|-8118 is of similar quality and for it we also 
only use scales Av > 32 km/s. 
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4. Hydrodynamical simulations 



4.1. 



The simulation 



For the hydrodynamic simulations, we use a Eulerian code based on a weighted essentially 
non-oscillatory (WENO) sciieme for a iiyperbolic system (Jiang & Shu, 1996; Shu, 1998). WENO 
realizes the idea of adaptive stencils in the reconstruction procedure based on the local smoothness 
of the numerical solution to automatically achieve high order accuracy and non-oscillatory 
property near discontinuities. It is extremely robust and stable for solutions containing strong 
shock and complex solution structures. This is especially important in studying intermittency. In 
the context of cosmological applications, we developed a hybrid N-body/hydrodynamical code 
that incorporates a Lagrangian particle-mesh algorithm to evolve the collision-less matter with 
the fifth order WENO scheme to solve the equation of gas dynamics (Feng et.al, 2002). 

In the present application, we computed the cosmic evolution of the coupled system of both 
dark matter and baryonic matter in a flat low density CDM model (ACDM), which is specified 
by the density parameter = 0.3, cosmological constant = 0.7, Hubble constant h = 0.7, 
and mass fluctuation within a sphere of radius 8h~^Mpc, ag = 0.9. The baryon fraction is fixed 
using the constraint from primordial nucleosynthesis as rij, = 0.0125h~^ (Walker et.al, 1991). The 
linear power spectrum is taken from the fitting formulae presented by Eisenstein &; Hu (1998). 
We perform the simulations in a 128'^ grid with an equal number of dark matter particles. In 
order to examine the effects of numerical resolution, two simulations on a periodic, cubical box of 
size 6h~^Mpc and 12h^^Mpc are generated. The simulations start at a redshift z = A9 and the 
results were outputed at the average redshifts of the each observation samples (for q0014-|-8118, 
Zm = 2.915, and HS1700-I-6416, Zm = 2.329). The timestep is chosen by the minimum value among 
three time scales. The first is from the Courant condition given by 



where Ax is the cell size, is the local sound speed, Ux, Uy and Uz are the local fluid velocities 

and cfl is the Courant number, typically, we take cfl = 0.6. The second constraint is imposed by 
cosmic expansion which requires that Aa/a < 0.02 within single time step. The last constraint 
comes from the requirement that a particle move no more than a fixed fraction of cell size in one 
time step. 

Atomic processes including ionization, radiative cooling and heating are modeled as in Cen 
(1996) in a primeval plasma of hydrogen and helium of composition (X = 0.76, Y = 0.24). The 
uniform UV-background of ionizing photons is assumed to have a power-law spectrum of the form 
J{i') = J21 X 10~^-'^(i'/i'if/)~"ergs~^cm~^sr~^Hz~-'^, where the photoionizing flux is normalized by 
parameter J21 at the Lyman limit frequency vhi, and is suddenly switched on at 2; ~ 6 to heat the 
gas and reionize the universe. 



6t < 



cfl X a{t)Ax 



(24) 



max(|uj;| -I- Cs, \uy\ + Cg, \uz\ + Cg) 



One-dimensional fields are extracted along randomly selecting lines of sight in the simulation 
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box. The density, temperature and velocity of the neutral gas fraction on grids are firstly Gaussian 
smoothed using FFT techniques which form the fundamental data set. The one-dimensional 
grid containing the physical quantities is further interpolated by a cubic spline. Using this 
one-dimensional grid, the optical depth r at each pixel is then obtained by integrating in real space 
and we include the effect of the peculiar velocity and convolve with Voigt thermal broadening. 
To have a fair comparison with observed spectra, r was Gaussian smoothed to match with the 
spectral resolutions of observation. The transmitted flux F = exp(— r) is normalized such that the 
mean flux decrement in the spectra match with observations. 

Each mock spectrum is sampled on a 2^° grid with the same spectral resolution as the 
observation. As the corresponding comoving scale for 2^^^ pixels is larger than the simulation box 
size, we replicate periodically the sample. To achieve the greatest the statistical independence, 
we randomly change the direction of line of sight while crossing the boundary of the simulation 
box. For each observed spectrum, 1000 mock spectra are generated. In addition, we also 
output 100 samples of the density distribution of both dark matter and neutral hydrogen on the 
one-dimensional grid. 

4.2. Power spectrum of hydrodynamical simulation samples 

Figure 2 compares the DWT power spectra measured in the mock samples and the 

observed data of HS1700. The comoving wavelength k in Fig. 2 is given by fc = 2tt/L, where 
L = (At;/i7o)[f^m(l + ^m)^ + ^a]~^^^- The 1-a errors for the simulation samples are derived from 
the 1000 mock spectra. For the observed spectrum of HS1700, we divide the first 2^"^ pixels of 
the observed spectrum into 16 subsections, each with 2^° pixels. The DWT power spectra are 
obtained for each subsection respectively, and the 1-a error is estimated from the variance of 
the 16 sections. [For an intermittent field, the error estimation of power spcctrTim is sensitively 
dependent on how one treats the intermittency (Jamkhcdkar, 2002). We do not discuss the details 
of this point here as it will not affect the following discussion.] 

Figure 2 shows that simulation samples with the 12 Mpc box are in good agreement with 
observations on all scales A; < 15 h Mpc^^, which corresponds to j = 10 in Fig. 1. The simulation 
samples with the 6 h~^Mpc box have higher power than observed sample on small scales, but 
lower power on large scales. This is expected because the simulation in the 6 h~^Mpc box lacks 
perturbations on scales A; < 1 h Mpc~^ and the normalization of the mean flux decrement will 
lead to the higher power on small scales to compensate for the loss of powers on large scales. 

On small scales /c > 15 h Mpc~^, the simulated power is significantly lower than the 
observation data. This shows again that noise contamination becomes important on scales A; > 15 
h Mpc~^, or j > 10. Therefore, as in §3.3, we use the simulation data only on scales A; > 15 h 
Mpc"\ or j < 10, or Av > 32 km/s. 
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5. Statistical discrepancy between baryonic gas and dark matter 

5.1. Relation between IGM and dark matter mass fields 

It is well known that in the linear approximation the density contrast of the IGM, (5{,(x, 
is equal to that of dark matter 6DMi^,t) on scales larger than the Jeans length. This is true 
even if initially ^^(x, tg) 7^ <^£)m(X)^o) (Bi, Borner & Chu 1993, Nusser 2000, Matarrese & 
Mohayee 2002). In other words, the density fluctuations of the IGM trace the dark matter fields 
point-by-point. Since the clustering of the IGM is governed by the gravity of dark matter, the 
IGM can be considered as "passive substance" with respect to the underlying dark matter. Thus, 
the point-by-point relations p5(x, t) oc p£)jv/(x, t) and 5b{x,t) = 5DM{^,t) on scales larger than 
the Jeans length also looks reasonable in the nonlinear regime. In this case, all the statistical 
properties of the IGM field on scales larger than the Jeans length are the same as the underlying 
dark matter field. 

However, it is known that a passive substance in a hydrodynamic system does not always 
behave like the underlying dominant mass field when non-linear evolution takes place. The 
statistical properties of a passive substance can decouple from those of the underlying field in 
the nonlinear regime. Moreover, the evolution of the IGM driven by the gravitational field of 
dark matter is described by a non-linear equation which is essentially the same as that used 
for modeling the statistical decoupling between a passive substance and its underlying field (see 
Appendix A.) Hence, we expect to see a statistical discrepancy between ^^(x, and 6DM{^,t) on 
scales larger than the Jeans length. 

To determine whether any statistical decoupling has occurred, we begin by assuming that 
5b(x, = (5dm(x, t). Then wavelet transforming yields e^^ = ej^^, and therefore 

(^i/^~r) = ([^1^]') = ([^ri')- (25) 

([e^^]^) and ([ej^^]^) are, respectively, the DWT power spectrum of the IGM and dark matter 
fields. Therefore, the first conclusion of the point-by-point relation (5f,(x, t) = SuM{^,t) is that the 
power spectra of the IGM and dark matter should be the same. Figure 3 indeed shows that the 
power spectra of the IGM and dark matter fields are about the same. 

With eq.(24), it is obvious that the power spectrum of the field (^^(x, t) — (5dm(x, t) should be 
equal to zero, or considering noise, very small. Figure 4 gives the power spectra ([e^; — ej;*^]^) and 
([e^J^) and shows that on large scales j < 8 {Av > 128 km/s) ([e^^ — ej^^]^) is much less than 
([e^J^). The ratio of ([e^j — ej^-^]^) to ([e^;]^) is < 0.05, which probably is due to noise. However, 
on scales j > 8 {Av > 128 km/s), this ratio gradually increases to ~ 0.4. This difference between 
large and small scales cannot be explained by the noise, as can be seen from Fig. 2 in which the 
error bars in the range 1 < /c < 10 h^^ Mpc (corresponding to j = 7 to 10) are about the same. 
Figure 4 indicates that (ej^ej]^^) < ([Cj;]^) on small scales, i.e., the point-by-point correlation 
between the IGM density fields and dark matter fields is weaker on small scales. This is a sign of 
the statistical discrepancy between the IGM and dark matter mass fields. 
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We can further make the case by calculating the structure functions for the field 

Sh{x,t) — (5£)m(x, t). The result is plotted in Fig. 5 and shows that for j = 6 the ratio of the 
moments is Gaussian, i.e. In S'|"/(5'|)" = ln(2n - 1)!! [eq.(5)]. On scales j > 8, In S'|"/(S'|)" 
is significantly larger than a Gaussian distributions which strongly indicates of the statistical 
discrepancy between 5b(x, t) and (5dm (x,t) is not due to noise. For a given n, In S'|"/(S'|)" 
increases with scale j so that the intermittent exponent ( [eq.(21)] is non-zero and negative, and 
therefore, the field (5fe(x, t) — 5dm{'^-, t) is intermittent. This result is consistent with the dynamical 
equation (A14), which requires that the difference (5;,(x, t) — 5dm(x, t) generally is an intermittent 
field. Thus, we conclude again that there is statistical discrepancy from the point-by-point relation 
6b{yi,t) = 5DM{^-,t). 

5.2. Relation between HI and dark matter 

The statistical discrepancy in the relation between the density fields of neutral hydrogen 
and dark matter is important as well. For instance, the pseudo-hydro technique assumes that 
dHi{y^,t) (X 5%M{yL,t) with 1.5 < a < 1.9 (Hui & Gnedin, 1997). This relation yields (ejf^ej^) > 0, 
and therefore the power spectrum ([ej^'^ + (fi^^Y) should be larger than — ^fi^Y)- Figure 6 
presents the power spectra {[e^^^ ± ^fi'^]'^) and shows ([Cjf^ + ej;^]^) > {[^fi^ — ^fi^Y) on large 
scales, but this difference is significantly smaller on smaller scales. 

Figure 7 gives the power spectra {[fji^ — f-fi^^) and ([e^J^) and reveals that it is similar to 
Fig. 4. All the results indicate that the positive correlation {^fi^^ji^^) > becomes weaker on 
scales J > 8 or A; < 3 h Mpc~^ or L > 2 h~^ Mpc. Since intermittent features become significant 
on scales j > 8, the dynamical assumption that pHi{^) oc p|)^(x) is applicable only for lower 
order statistics (like the power spectrum.) 

6. Comparison of intermittent features between real and simulation samples 

6.1. Structure functions of real and simulation samples 

We now calculate the structure functions for the transmitted flux fluctuations of the simulation 
samples. The results are illustrated in Fig. 8 for HS1700 and Fig. 9 for q0014. The error bars 
are the 1-a deviations estimated from the 1000 simulation samples. Comparing Figs. 8 and 9 
with Fig. 2, we see that lnS'|"/(S'|)^ with n = 2 has smaller errors than the power spectrum. 
Even at 6th and 8th orders, the j-dependencc of In S'|"/(S'|)"' can be very well fitted by a line. 
This implies that to measure the higher order statistical behavior of an intermittent field, the 
structure functions arc effective and stable. For an intermittent field, the large uncertainty in the 
power spectrum are caused by rare and improbable high spikes in the fluctuations. However, the 
structure function is an ensemble average of the ratio between S?" and (S?)" and this reduces the 
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effect of individual high spikes. 

Figures 8 and 9 also show that for a given n, In Sj"'/{Sj)^ increases with scale j so that the 
intermittent exponent ^ [eq.(21)] is non-zero and negative. The transmitted flux fluctuations axe, 
hence, intermittent. In Fig. 10, we show In 5|"'/(5'|)"' vs. j of HS1700 mock samples produced by 
both the 12 and 6 h~-'^Mpc box. The 12 h^^Mpc samples is a little higher In S'?"/(S'|)"', but there 
arc no significant differences between the simulation samples on scales considered. 

The structure functions measured for the real data of HS1700 is also shown in Fig. 10, in 
which the error bars HS1700 are found by bootstrap re-sampling. Within the error bars, the 
observed structure function is consistent with the 12 Mpc box simulation samples for the LCDM 
model. Although the consistency shown in Fig. 10 is still not perfect, it is substantially better 
than the results obtained from the pseudo hydro simulations. Pando et al (2002) found that the 
structure functions of the pseudo hydro simulation samples were generally larger than those of the 
real data [see Fig. 7 of Pando et al. (2002)]. 

The structure functions for both HS1700 and q0014 have similar behavior although they lie 
at different redshifts. The result demonstrates that one cannot see a significant redshift evolution 
in the range of 2 < z < 3 with the j-dependence of the structure functions. 

6.2. n-dependence of the structure function 

We now turn to the n-depcndcncc of the structure functions. Figures 11 and 12 are, 
respectively, ln5'?"/(5?)"' vs. n for the simulation samples and real data of HS1700 and q0014 
on scales j = 9 and 10. The figures show again that the simulation result is consistent with 
observations on the scales considered. 

For a Gaussian field, the n-dependence of structure function is lnS'|"/(5|)" = ln(2n — 1)!! 
[eq.(5)], which is also plotted in Figs. 11 and 12. We note, from Figs. 11 and 12, that although 

hiS'|"/(S'|)" is above that of a Gaussian field, In Sj"- / (Sj)"- seems to "parallel" the Gaussian curve. 
This means that although the field on scale j = 9 is non-Gaussian, the PDF of ej*^ doesn't have a 
significantly long tail. On scale j = 10, the slope of S'|"/(S'|)" becomes larger than (2n — 1)!!, 
which indicates that the PDF of ej^ tends to be long tailed. 

The n-dependence of ln5'?"/(5|)" appear to be more sensitive to redshift evolution. Figure 

11 shows that the relation of 5|"/(5'|)" for j = 9 is clearly above the Gaussian field, while for Fig. 

12 it is only marginally different from a Gaussian field. Thus, the intermittent behavior on scale 
j = 9 probably developed at redshift z < 3. 

Similar to the power spectrum (Fig. 1), the results of Figs. 8-12 result is /-independent on 
the entire scale range considered for / = 1 to 3 (Pando et al 2002). To see the effect of metal lines, 
we plot Fig. 13, which shows the n-dependence of In5'?"/(S'|)"' {j = 10) with and without metal 
line removal. The result for j = 10 is less affected by the metal line removing because most metal 



lines have width less than Av = 32 km/s {j = 10). 



6.3. Intermittent exponent 

The n-dependence of the intermittent exponent ( can be calculated by eq.(17), which gives 

1 S'^"' 

C(n) = -^ln2— |-+const, (26) 
J [^j ) 

or 

1 

CN-C(l) = -^ln2^. (27) 

For the simulation samples, the relation of ln2 S'j'^ / {S'j)'^ vs. n shown in Fig. 10 can be very well 
fitted by a line in the range from j = 7 to 10 (or 2 < A; < 16 h Mpc~^), and therefore, one can find 
a constant C{n) for each given n. The real data is not as well fitted over the same scales, but we 
nonetheless show the results in figures 14 and 15. 

Again the full hydro simulation shows better result than the pseudo hydro simulation. The 
value of I CI for the real data is found to be systematically lower than that given by pseudo hydro 
simulation for the LCDM models (Pando at al. 2002). Although, the value \C,\ from the hydro 
simulation is still higher than real data, the differences are no longer larger than the error bars for 
all n. Figure 16 shows the intermittent exponent for samples with 12 and 6 h~^Mpc boxes. Note 
that the conclusions drawn from Figs. 14 and 15 hold also for samples with the 6 h~^Mpc box. 



6.4. Intermittency of HI density field 

It has been argued that the intermittency of the Lya transmitted flux field F{\) may not 
indicate intermittency in the IGM distribution (^^(A) or HI distribution 6hi{^)- The reason is that 
in the exponential relation between the flux and optical depth -F(A) = Fce^'^^'^\ and the optical 
depth r(A) is approximately related to ^^(A) by 

t(A) = A\phiX)/pbr = A[l + 6hiX)r, (28) 

where the parameter A depends on the cosmic baryonic density, the photoionization rate of HI, 
and the temperature of IGM. Therefore, even when Si,{X) or 5hi{X) are Gaussian, the PDF of the 
flux -F(A) might be lognormal due to the exponential relation between F{X) and 6hi{X). 

However as has been emphasized in §2, the intermittency of F{X) field is not measured by the 
PDF of F(A), but its difference AsxF{X) = F{X + 5X) - F{X). The structure functions are defined 
by the density difference. Thus, we have approximately 



AsxFiX) ~ -e-^WAar(A) 



(29) 
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where Asxt{X) = r(A + SX) — r(A). Using eq.(27), we have 

Aar(A) ~ -Aa[l + 6Hi\)r~'Asx5b{\), (30) 

where Ag^SbiX) = Sb{X + SX) - 6b{X). 

Eqs.(29) and (30) show that if the transmitted flux fluctuations are intermittent, the baryonic 

matter density field and HI field should also be intermittent (Zhan, Jamkhedkar, &; Fang 2001). 
Figure 17 presents the structure functions for the HI field of the 12 Mpc box samples simulation 
and shows that the HI field is indeed intermittent. In fact, this field has larger structure functions 
than the flux fluctuations probably because the transmitted flux cannot trace the intermittency of 
optical depth at heavy absorption regions. 



7. Discussion and conclusions 

We analyzed the intermittent behavior of the Lya transmitted flux fluctuations of real data 
and samples produced by full hydro simulations. This analysis covers scales from ~ 10 h""^ Mpc 
to 400 h~^ kpc (j = 6 to 10), and statistical orders 2n from 2 to 8. The intermittent behavior is 
found to be significant on scales less than about '-^ 1 h^"*^ Mpc, and rcdshift z < 3. Over the entire 
scale range and statistical orders considered, the intermittent behavior of the real data matches 
the samples produced by full hydrodynamic simulations in the LCDM model. 

We have shown that the point-by-point correlation of the IGM and HI fields with dark matter 
field becomes weaker on scales less than about 2 h~^ Mpc, which is much larger than the Jeans 
length at the mean density. Although the power spectrum of the IGM field is still about the same 
as that of dark matter at these scales, the mass density field of the IGM is no longer point-by-point 
proportional to the underlying dark matter field. 

By using intermittency we have detected a the dynamical discrepancy between the IGM 
field and the dark matter distribution. The study of cosmic large scale structure has mainly 

concentrated on either the two-point correlation function and its scaling, or on the massive 
collapsed halos of dark matter. Intermittency may bridge the gap between the two approaches 
and provide new, physically interesting insights in cosmic clustering. 
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A. Dynamical equations of the IGM and dark matter mass field 

A.l. Basic equations 

Let us consider a flat universe having cosmic factor a{t) oc t^/^ and dominated by dark 
matter. We describe tlie dark matter by a mass density field p{x, t) and a peculiar velocity field 
v(x, t), where x is the comoving coordinate. In hydro dynamical descriptions, the equations of 
dark matter consist of the continuity, the momentum, and the gravitational potential equations as 
follows (Wasserman 1978) 

^ + -V-(l + (5)v = (Al) 
at a 

^ + (V ■ V)v = - (A2) 

= ^T^Ga^pS (A3) 

where density perturbation (5(x, t) = [p(x, t) — p]/p, and the mean density p = l/GirGt^ oc a"^. 
The peculiar gravitational potential cf) is zero (or constant) when density perturbation 5 = 0. The 
operator V is acting on the comoving coordinate x. For the growth mode of the perturbations, 
the velocity is irrotational. In this case, one can define a velocity potential by 

V = --Vip (A4) 

a 

With this potential, the momentum equation (A2) can be rewritten as the Bernoulli equation 

t-^(V.f = *. (A5) 

Since there is only gravitational interaction between dark matter and the IGM or cosmic 
baryonic gas, it is convenient to describe the IGM by its mass density field P6(x, t) and velocity 
field Vb(x, t). The hydrodynamical equations of the IGM are 

^ + lv-il + d,)v, = (A6) 

^ + (V5 • V)v, = --Vpb - V<P (A7) 
dt pb 

where the density perturbation of the IGM 55(x, t) = [pi,{x,t) — pbj/pb, and pb the mean density 
of the IGM. The gravity of the IGM is negligible. The evolution of the IGM is governed by the 
gravitation of dark matter, and therefore, the gravitational potential (f) in eq.(A7) is still given by 
eq. (A3). 

To sketch the gravitational clustering of the IGM, we will not consider the details of heating 
and cooling. Thermal processes are generally local, and therefore, it is reasonable to describe the 
thermal processes by a polytropic relation pi,{'x.,t) oc p2'(x, t). Thus eq.(A7) becomes 

^^ + (v6-V)v6 = (A8) 

ot prup (1 + Ob) 



-19- 



where the parameter /x is the mean molecular weight of IGM particles, and nip the proton mass. 
In this case, wc don't need the energy equation and the IGM temperature evolves as T oc p''~^, or 
T = To(l + 5h)^^^. Eq.(A8) is different from eq.(A2) only by the term with temperature T. If we 
treat this term in the linear approximation, we have 

where ipb is the velocity potential for IGM field defined by 

V6 = --Vifb. (AlO) 
a 



The coefficient is given by 



where D{t) describes the linear growth behavior. 

Eq. ( A9) shows that the clustering of the IGM field is completely controlled by the gravity of 
underlying mass fields of dark matter. It can then be considered as a "passive substance" with 
respect to the underlying dark matter. The term with ui, in eq.(A9) acts like a viscosity which is 
due to thermal diffusion characterized by the Jeans length kj = {a^ /t'^){ump/^kBTo). In linear 
regime, eqs.(A5) and (A9) yield solution ipf, = (p on scales larger than the Jeans length. 



A. 2. Equation of the difTerence between the IGM and dark matter fields 
To study the possible difference between the IGM and dark matter, we define a variable 9 by 

e = ipb-^ (A12) 

which describes the deviation of IGM from dark matter. With eqs. (A5) and (A9), we have 

When = 0, this has solution = const, if the initial condition is = const. During linear 
evolution, we have = 0, and therefore = is the initial condition for eq.(A13). Thus, the 
solution 6 = is correct even in the nonlinear regime if z/;, = 0. There is no deviation of IGM 
from dark matter. 

In the case of f;, 7^ 0, = will no longer be the solution of eq.(A13) even if it is = 
initially. Considering that the potential 93 of dark matter essentially is stochastic, the r.h.s. of 
eq.(13) plays the role of stochastic forces. It drives to be non-zero. This leads to a statistical 
or stochastic discrepancy of the IGM field from dark matter field. This discrepancy is not trivial 
even when we consider the linear approximation of eq.(A13) with respect to 0, i.e. 

^_^,ve,(v.)-5v=e = 5vV (am, 
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Eq.(14) essentially is the same as that widely used to dynamically model the discrepancy, or 
decoupling, of the statistical properties of the passive substance from the underlying field (for a 
review, Shraiman & Siggia, 2001). Since the advection term ^(VG)(V(^) doesn't depend on the 
diffusion scales, the statistical discrepancy of the passive substance from the underlying field can 
appear on scales larger than the diffusion scale. 

The field O given by eqs.(A13) and (A14) is generally intermittent. It has been shown that 
the field @ given by eq.(Al4) will be intermittent, even when the field if is Gaussian (Kraichnan, 
1994). Eq.(A13) is similar to the stochastic- force-driven Burgers' equation, or the so-called KPZ 
equation (Kardar, Parisi, & Zhang 1986, Berera Sz Fang 1994, Barabasi, &; Stanley, 1995, Jones 
1999), which is typical of the dynamical model of intermittency (Polyakov 1995; Balkovsky et al. 
1997; E et al. 1997). Therefore, the difference between the IGM and dark matter fields, 9, is 
probably intermittent. 
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Fig. 1. — The power spectra of the HS 1700+6416 transmitted flux for 1.) the fitted flux (triangle), 
2.) conditional counting with / = 1 (square), and 3.) with / = 3 (pantegan). The physical scale 
for j is Av = 2^^~^ x 4 km/s. 




Fig. 2. — The power spectra of HS1700+6416 for 1.) real data (square); 2.) hydro simulation in a 
12 Mpc box (circle), and 3.) in 6 Mpc box (triangle). The comoving wavenumber k is calculated 
in the LCDM model. For clarity, the power spectra of the simulations are shifted slightly. 
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Fig. 3. — Power spectra of IGM mass field 6f, (square) and dark matter Sdm (circle) from the hydro 
simulation in the 12 Mpc box. The mean and variance are calculated from the 100 realizations. 
The physical scale for j is Av = 2^^~^ x 4 km/s. For clarity, the power spectrum of the IGM is 
shifted slightly. 
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Fig. 4. — Power spectra of 5b (circle) and 5b — 5dm (square) of the hydro simulation samples for 12 
Mpc box. The physical scales for j is Av = 2^^~^ x 4 km/s. 
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Fig. 5. — Structure functions ln2[S'|"/ (Sj)] vs. j of the difference Sb — Sdm of the hydro simulation 
samples in the 12 Mpc box. The physical scale for j is At; = 2^^~^ x 4 km/s. 
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Fig. 6. — Power spectra of 6hi + Sdm (circle) and 6hi 
simulation samples. The physical scale for j is Av = 
spectrum of 6hi — Sdm is shifted slightly. 



— Sdm (square) of the 12 Mpc box hydro 
^ 4 km/s. For clarity, the power 
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Fig. 7. — Power spectra of 5hi (circle) and 5hi — Sdm (square) of the 12 Mpc box hydro simulation 
samples. The physical scales for j is Av = 2^^~^ x 4 km/s. For clarity, the power spectrum of 
Shi — Sdm is shifted slightly. 
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Fig. 8. — Structure functions ln2[S'|"/ (5'|)] vs. j of the hydro simulation samples for HS1700+6416 
in the 12 Mpc box. The physical scale for j is Av = 2^^~^ x 4 km/s. 
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Fig. 10.— Structure function ln2[S]"- / {S])] vs. j for HS1700+6416 of real data (circle) and hydro 
simulation with 12 (triangle) and 6 (square) Mpc box. The physical scale for j is Av = 2^^^^ x 4 
km/s. For clarity, the results of the two simulation samples are shifted slightly. 
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Fig. 11. — In vs. n for simulation sample (circle) and real data (square) of HS1700+6416, 

and j = 9 and 10, or Av = 64 and 32 km/s. The dash line is for a Gaussian field. For clarify, the 
results of HS1700 data are shifted a little slightly. 
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Fig. 13. — The same as Fig. 11 for j = 10. The top panel is for real data with metal line removing. 
The bottom panel is without metal line removing, which has about the same error bars as the top 
panel. The dash line is for a Gaussian field. For clarify, the results of HS1700 data are shifted a 
little slightly. 
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Fig. 14. — Intermittent exponent C vs. n for simulation data (circle) and real sample (square) of 
HS1700. For clarity, the results of HS1700 data are shifted a little bit to right. 
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Fig. 15.— 



The same as Fig. 11, but for q0014+8118 
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Fig. 16. — Intermittent exponent ^ vs. n of HS1700 samples given by hydro simulation in 12 
(square) and 6 (circle) Mpc box. For clarify, the results for 12 Mpc box are shifted a little to the 
right. 
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Fig. 17. — Structure functions ln2[iS'|'^/(-S'?)] vs. j of HI density field Shi- The physical scale for j 
is Av = 2^^~^ X 4 km/s. 



